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source code implicit finite element solver 



rate/temp/press dependent von mises isotropic plasticity 

umat for abaqus 5-5- nonlinear strain hardening- 

2d/3d problems with the exception of plane stress 

by omar a hasan last modified 05-02-^ 

user must specify differential hardening data in umat 

and dimension hardening table appropriately 

must have atleast two sets of points in table 

subroutine umat C stress Tstateviddsddeisseisptiiscd-i 

1 rpl nddsddt-,drplde-,drpldt ^ 

2 stranidstranitimeTdtimeitempidtempipredef idpredicmnamei 

3 ndi ^nshr-,ntens -ins tat Vi props nnprops i coords tdrot-i pnewdti 

4 celent idfgrdQ n dfgrdl-.no el -i npti layer Tksptnkstep-r kinc) 

include ' aba_param- inc ' 
character*fl cmname 

dimension stross(ntans),5hf3tev(nstatv) i 

1 ddsdde (ntens ntens) t ddsdd t ( ntens ) -i dr p Ide ( nt ens ) -< 

2 stran (ntens) dst ran ( ntens ) n time (2) -.predef (1) -.pred(l) n 

3 props (npr ops) -.coords (3) -.drot (3-iH) i df grdO ( 3 3 ) -> d f gr dl (3 -> 3 ) 

dimension flow(L) 

parameter (zero = D-dQTone=l.dDn two=2.dDithree = 3- do-. si x=t-d0n 
1 newton = tDi toler = l .Qd-Sntwbth^D .bbbfctibbbbfcibdD) 



cannot be used for plane stress 



props(l) - e (Pa) (temperature dependent) 
props(2) - nu 

props(3) - rate sensitivity (temperature dependent) 
props(H) - intrinsic flow rate (temperature dependent) 
props(S) - prQssure sensitivity 

calls uhard for curve of intrinsic strength vs. plastic strain 



material properties 
emod=props (1) 
enu=props (2) 

ebulk3=emod/(one-two*enu) 

eg2=emod/(one+enu) 

eg = eg2/tu/o 

eg3=three*eg 

elam=(ebulk3-eg2) /three 

r lp2m= el am+eg2/ three 

r a tesf =pr ops ( 3 ) 

rrates=one/ratesf 

dtebsO=dtime*props(M) 
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source code implicit finite element solver 



psf =props (5) 

esi=dstran(l)**E+dstran(E) **E+ds tran C 3 ) **E 
do kl=nd i +1 t ntens 

esi=esi+two*(dstran(kl)/two)**E 
end do 

esi = sqrt(tujbth*esi) 

s_rate = max (l.d-lD-,esi/dtime) 

elastic stiffness 
call aset (ddsddeTzerointens*ntens) 
do kl = l-fndi 
do kE = l-.ndi 

ddsddeCkEnkl)=elam 
end do 

ddsdde(klnkl)=egE + elani 
end do 

do kl = ndi+l-jntens 

ddsdde(klnkl)=eg 
end do 

recover equivalent plastic strain & equivalent stress 
and hydrostatic stress at start of step 

eqplas=statev CI) 

qold=statev(E) 

hydr_o= (stress(l)+stress(E)+stress(3))/three 

calculate predictor stress 
do kl = l-,ntens 
do kE = l->ntens 

stress (kE) =s tress (kE) +ddsdde (kS-ikl) *dstr an (kl) 
end do 
end do 

calculate equivalent von mises stress 

smises=Cstress(l)-stressCE))**E+Cstress(E)-stress(3))**E 
1 + (stress (3) -stress (1) )**E 

do kl=ndi+ln ntens 

smises=smises+six*stress(kl) **2 
end do 

smises=sqrt(smises/two) 

get differential hardening from the specified hardening curve 
call uhard (syielDThardieqplas) 

determine if actively yielding 
if (timed) -gt-D-dD) then 

separate the hydrostatic from the deviatoric stress 
calculate the flow direction 
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source code implicit finite element solver 



shydro=(stress(l)+5tress(E)+stress(3))/three 
do kl = l->ndi 

flow(kl)=(stress(kl)-shydro)/smises 
end do 

do kl=ndi+l intens 

flow(kl)=stress(kl)/smises 
end do 

c 

c solve for equivalent von ntises stress 

c and equivalent plastic strain increment using newton iterati 

on 

syie ld = sy i elO 

c use this to minimize iterations during elastic deformation ( 

1) 

c deqpl=dtebsO*exp((smises-syield)*ratesf) 

c use this to minimize iterations during plastic deformation ( 

E) 

deqpl=esi 

do kewton=l t newton 

deqpl=max(deqplnl-d-5D) 

qhs=smises-eg3*deqpl-syield-rrates*dlog(deqpl/dtebsO) 
rhs=qhs+psf*shydro 

deqpl=deqpl + deqpl*rhs/(deqpl*(eg3+hard)+rrates) 
call uhard(syieldihardTeqplas+deqpl) 
if (abs(rhs) - It . toler*bQ - dO) goto 10 
end do 

write(7iE) newton 
E format (// n 3E3x i 1 ***warni ng - plasticity algorithm did not 

1 'converge after 'ii3i f iterations') 

write(7-i*)dstran(l)-idstran(E)-idstran(3)-.dstran(M) 
write (7-i*0 dstran (5) Tdstran CD TesiismisesnStatev(l) 
write (?-,*) statev(E) -,statev(3) istatev(4) -,statev<5) 
write(7n*)qhsideqplirhsnshydronStress(l)TStress(E) 
write(7-i*Ostress(3)TStress(H)nstress<h)iStress(b) 
ID continue 

c 

c the new equivalent deviatoric stress (q) is 

q=sy ield+rrates*dlog (deqpl/dtebsD) -psf *shydro 

c 

c update stressi elastic and plastic strains and 

c equivalent plastic strain 

do kl=l-«ndi 

stress(kl)-flow(kl)*q+shydro 

end do 

do kl=ndi+l intens 

stress (kl)=flow(kl)*q 
end do 

eqplas=eqplas+deqpl 
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source code implicit finite element solver 



c calculate plastic dissipation 

spd=deqpl*(qold+q)/two 

c 

c formulate the jacobian (material tangent) 

c first calculate effective moduli 

effg=eg*q/smises 
e f f gE=two*ef f g 
effg3=three/two*effgE 
efflam=(ebulk3-effgE) /three 
hardl=hard+rrates/deqpl 
effhrd=eg3*hardl/(eg3+hardl)-effg3 
cee=-ebulk3*psf*eg*deqpl/smises 
do kl = l-indi 
do kE=lindi 

ddsdde(kEnkl)=ef f lam+cee*f low(kE) 
end do 

ddsdde (klikl)-ef f gE + ef f lam+cee*f low Ckl) 
end do 

do kl=ndi+lintens 

ddsdde(kl-,kl)=ef f g 
end do 

do kl=lintens 
do kE=lintens 

ddsdde (kS-ikl) =ddsdde ( kE n kl) +e f f hrd*f low (kE) *f low (kl) 
end do 
end do 
endi f 

c 

c store state variables in array 

c equiv s tr ai n -i mises stress *• plastic strain ratenelastic strain 
c rate and iterations to convergence 

statev(l)=eqplas 

state v CE ) =q 

statev(3)=deqpl/dtime 

statev(M)=esi/dtimQ 

statev(5)=kewton 

c 

return 
end 

-c 

subroutine uhardCsyieldihardneqplas) 

c 

include ' aba_param . i nc ' 
c table must be dimensioned correctly below: 

dimension table(Ei7) 
parameter(zero=0-dO) 
c nbv 313 hardening table 

nva lue=7 
c this is room temp data 

table(lnl)=0-DQda 



/VP 
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source code implicit finite element solver 



table(2-il)=0-0 

table(l-.2)=-5.5 c !5d0 

table(S-.S)=D-151 

table(l-.3)=-3.D4dO 

table(5-,3)=D.337 

tabled. M) = M.7ELdO 

table(eiM)=0.SM5 

table(li5)*lM-mdO 

table(2i5)=Q.73b 

table (lib) =Mfi -IMbdD 

table(2-ife)=3i.CH3 

table(ln7)=270i4-4dD 

table(a,7)=17.0fib 

do kl = l -tnvalue-1 

eqpll = table(2-.kl + l> 
if(eqplas.lt-eqpll) then 
eqplD=table(2ikl) 
current yield stress and hardening 

deqpl=eqpll-eqplD 
syielD-table (likl) 
syiell = table (1-ikl+l) 
dsyiel=syiell-syielO 
hard^dsyiel/deqpl 

syield = syielQ + '(eqplas-eqplD)*hard 
goto ID 
endi f 
end do 
10 continue 

return 
end 
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source code explicit finite element solver 



c vectorized user material subroutine for shell and plane 

c stress elements (abaqusS-5) 

c rate/temp dependent isotropic plasticity with linear 

c elasticity! strain softening/hardening & press- depnd- 

c yield 

c by omar a hasan (hasanScrd-ge.com) 

c last modified 05-03-^L 



c 



c 

ps) 



subroutine vumatC 

read only variables C unmodi f i abl e ) 

1 nblockindir-inshrinstatev-infieldvinprops-ilannealT 

2 step_time-» total_tin\e ndt -icmname ^coorOp ichar Jength^ 

3 propsidensity nstrain^inc t rel_spin_i nc i 

M temp_oldTStretch„oldidef grad_old-> f ield_oldn 

5 s tress_o 1 d -i state_ ol d t ener_intern_ol d i ener„inelas_old t 

L temp_newi stretch_new-i dG f grad_new-« f ield__newn 

write only variables (modifiable) 
7 stress_new-i state_newT ener_intern_new Tener„inel as„new) 

include 1 vaba^param - inc ' 

dimension coord_mp ( nblock i *) i char_length ( nblock ) -i props (npro 

1 density (nblock) n str ain_i nc C nb lock -> ndir + nshr ) t 

5 rel_sp inline (nblock-mshr) -i temp_old ( nblock ) n 

3 stretch_old £nblockindir+nshr)n 

M def grad_old (nblock-.ndir + nshr + nshr) if ield_old (nblocknnfieldv 

5 stress_old Cnblockindir+nshr) t state__old ( nblock t nstatev ) n 
b ener_i ntern_o Id ( nblock ) -i ene r_i ne 1 as_o Id C nblock ) -i 

7 temp_new( nblock) n str etch_new C nblock t ndir +nshr ) n 

6 def grad_new Cnblockindir + nshr + nshr) f i el d_ne w ( nb 1 ock n n f i e 1 d v 

^ stress_new ( nblock -mdir + nshr) i state_new ( nblock instatev)n 
1 ener_i nt ern_new ( nbl ock ) -i ener_inel as„new ( nbl ock ) 

integer limit 
parameter (limit=M0) 
dimension tableCBi^) 
character*fl cm name 

parameter(2ero=D-dDTone=l-dDTtwo=2-dDnthree=3*dDTsix=b-d0T 
1 f our = 4 -dDnoneptf = 1 - SdD -. zept = 0 - SSdD -i tubth = D - Lbt.bbbt.bbLbdD 
E eitee=aO-dD) 



c propsCl) - e- modulus (temperature dependent) 

c props(2) - nu- poisson ratio 

c 

c Properties 3 and M descibe the rate sensitivity of yield ba 

sed on a plot of 
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source code explicit finite element solver 



c yield stress (x-axis) vs ln(strain rate) y-axis 

c 

c propsC3) - rate sensitivity (temperature dependent) SLOPE 

c props(H) - intrinsic flow rate (temperature dependent) INTER 

CPT 
c 

c Property 5 descibes the pressure sensitivity of yield 

c 

c props(S) - pressure sensitivity factor 

c 

c Property b is the failure criterion either an equivalen 

t plastic strain 

c for ductile failure or a maximum principal stress for britt 

le failure 

c 

c props(b) - failure criterion 

c 

c NOTE -THESE FOLLOWING TldO LINES WOULD APPEAR IN THE ABA£2US 

EXPLICIT 



c INPUT DECK 

c 

c *USER MATERIAL -i C0NSTANTS = 5 

c 2 - EHeiiD - MO -i 3 . STe-?-.! . Mfle-IM -iD - lb 

c *DEPVAR-iDELETE = b 

c fl 

c 

c 

c material properties 



emod=props (1 ) 
enu=props (2) 

ebulk3=emod/(one-two*enu) 

eg2=emod/(one+enu) 

eg=eg2/two 

eg3=three*eg 

elam=(ebulk3-eg2) /three 

e lp2g=e lam+eg2 

ratesf =props (3) 

dtebsD=dt*props(M) 

ps f = props ( 5 ) 

rrates=one/ratesf 

f ailst=props (b) 

table(lil)=0.0 

table(2il)=D-D 

table(li2)=b-S 

table(2i2)=0.15 

tabie(ln3)=17-^3 

table(2i3)=0.3S 

table(l-,4)=3M.M7 

table(2->M)=0-5S 

table(l-iS) =53.DT 
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source cede explicit finite element solver 



table(BiS)=D.7S 
table (lib) =70.35 
table(5-,b)=G- c iS 
tabled-.?)^!. □! 
table(En7)=1.15 
tabled-,B)=14bdb 
table(Eifi)=l-35 
tabled-. T)=5Dl-3 
table(E-, c l)=1.5S 

c 

do 10D idinblock 

c 

c initialize state variables 

eqplas=state_old (id) 
s m_o 1 d = s t a t e_o 1 d ( i -i 2 ) 
icont=state_old ( i t 3) 
tstart = total__time-dt 
if ( tstart - It d-e-b) then 
icontd 

state_old ( i ^ L) = one 
endi f 

c 

if Cstate_old(i-,fc) -lt.D-5) then 
st ate_new (i-«b)=zero 
goto 100 
endi f 

c 

c get hardening modulus and intrinsic resistance at t 

hard=(tabled-iicont + l)-tabledticont))/ 
1 (table(2dcont + l)-table(2dcont)) 

s_intr = table d-»icont> +hard*(eqplas- table (2-iicont ) ) 

c 

c calculate predictor stress 

trace2 = str a inline ( i d ) +strain_inc ( i n 2 > 

del_e33=-elam*traceE/elp2g 

s igllo = stress_old (iil)+egB*strai n_i nc 

sigEEo=stress_old ( i tE ) +egE*s trai n_inc ( i n E ) 

sig33=zero 

sigl2 = stress__old C i t M ) +egE*str ain_inc ( i 1 4 ) 
ss!Es = six*(sigliE**E) 

c 

c since str ain_inc ( i 3 ) is not known apriori-i loop 3 

c times without checking for convergence (works very well 

c in practise by reducing sig33 to 0 - 00000Ql*syield ) 

do EDO iidi3 

trace=trace2+del_e33 

sigll=sigllo+elam*trace 

sig22=sig22o+elam*trace 

c 

c calculate equivalent von mises stress from deviatoric 
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source code explicit finite element solver 



c component of trial (predictor) stress- 

smises= (sigll-sigBE) **E+(sigBE) **E + (sigll) **E 

smises=smises+sslBs 

smises=sqrt(smises/two) 
c avoid division by zero during first iteration 

smises=max(oneismises) 

c 

c separate the hydrostatic from the deviatoric stress 

c calculate the flow direction 

shydro=(sigll+sigE2)/three 

flowll=(sigll-shydro)/smises 

flowE2=(sigEE-shydro)/srnises 

flou/33=(sig33-shydro)/smises 

flowlE=sigl2/smises 

c 

c solve for equivalent von mises stress and equivalent 

c plastic strain increment 

adfp=-psf*shydro*ratesf 

deqpl=dtebsO*exp(( sm_ol d-s_intr ) *ratesf +adf p ) 
s m n e w = s m i s e s - e g 3 * d e q p 1 

c 

c update e33 

opfe=oneptf*deqpl 
d_epll = opf e*f lowll 
d_epEE=opf e*f lowEE 
d_ep33 ss opf e*f low33 
d_eplE=opf e*f low IE 
d_eell=strain_inc C i 1 ) -d_epll 
d_eeEE =str a in_inc ( i E ) -d_ep EE 
d_ee33=-el am* ( d_ee ll+d_ee EE ) /elpEg 
d_eelE=strain_inc(i ) -d_eplE 
del„e33=d_ee33+d_ep33 
EOD continue 

esi=strain_inc ( i i 1) **E+str a inline ( i iE ) **E + 
1 del_e33**B+two*strain„inc C i iM ) **E 
esi=sqrt(esi*twbth) 
strain_inc C i t3) =del_e33 

c 

c update stressn equivalent plastic straini location 

c of plastic strain counter and state variables 

stress_new( i -il) = f 1 owll*sm_new+shy dr o 

stress_new ( i -» S ) = f low22*sm_new+shydro 

stress_new ( i 3 ) = zero 

stress_new ( i i M ) = f lowlB*sm_ne w 

eqplas=eqplas+deqpl 

if (eqplas*gt-table(Eiicont + D) icont = icont + l 
cstate_new( i 1 ) =s tate_ol d (i il) +d_eell 
cstate_new ( i ->2 ) = state_old ( i i E ) + d_eeES 
cstate_new(i 3 ) =s tate^ol d C i n 3 ) +d_eelE 
cstate_new < i -i M ) =state__ol d ( i i M) + d_epll 



& /£j? 
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source code explicit finite element solver 



cstate_new ( i t 5 ) =s tat e_old < i *> 5) +d_ep EE 

cs tate__new C i -rb ) = state_old < i ->b ) +d_eplE 
c save state variables: plastic strainn vm stress^ total 

c strain ratei plastic strain ratei failure criterion flag 

stat e — n ew(i-il)=eqplas 

state__new ( i t E ) =sm_new 

state_new (in3)=icont 

stat e_n ew(i-iM)=esi/dt 

state„new( i n5) =deqpl/dt 

state_new( iib) =state_old ( i i b ) 

c 

bee = - C str ess_new C i 1 1 ) +str ess„new ( i E ) ) 
beeE=bee*bee 

cee=stress_new ( i -i 1 ) *stress__new ( i t E ) -stress„new( i -i M ) * 
1 s tr e ss_n bw ( i •> 4 ) 

f root=bee2-four*cee 
f frot=max(oneifroot) 
sqbmMc=sqrt(ffrot) 
pmax=(-bee+sqbm4c)/two 
pmin=(-bee-sqbm l 4c)/two 
state_new (i->?)=pmax 
stat e_n ew(iTB)=pniin 
c UNPAINTED 

f ailst = & c i -Ob 

if (pmax.gt-failst) state__new(i ibS-zero 
strain based failure criterion 
if (eqplas-gt-failst) state — new (iib^zero 
c update plastic dissipation 

plastic_work_inc=deqpl * ( sm_old+sm_new) /two 
ener_inelas_new ( i ) -ener_ine 1 as_o Id < i ) + 
1 plastic__work_inc/density (i ) 

c 

1DD continue 
return 
end 



